Impact of complex boundary on the hydrodynamic properties of methane nanofluidic flow via non-equilibrium multiscale molecular dynamics simulation

Understanding the impact of complex boundary on the hydrodynamic properties of methane nanofluidic is significant for production optimization and design of energy-saving emission reduction devices. In the molecule scale, however, the microscopic mechanisms of the influence of the complex boundary on the hydrodynamic characteristics are still not well understood. In this study, a mixture boundary Poiseuille flow model is proposed to study the hydrodynamic properties and explore the molecular mechanisms of confined methane nanofluidic using the Non-equilibrium multiscale molecular dynamics simulation (NEMSMD). In order to investigate the influences of nonslip and rough boundary on hydrodynamic behavior of nanofluidic by the present model in one simulation, the coordinate transformation methods regarding the local symmetry is showed. Simulation results show that the atom number density, velocity and temperature profiles present significant differences near the nonslip boundary and rough wall surface. Moreover, the slip length of methane nanofluidic near the rough boundary decreases with the increasing of the temperature. Furthermore, the viscosity values are calculated by parabolic fit of the local velocity data based on the present model, which demonstrates that the impact of the nonslip boundary on the shear viscosity compared with the experiment result is less than one obtained using the rough boundary. In addition, the local contours of rotational and translational energy are plotted, which show that the rotational and translational energies of nonslip boundary are obvious higher than those of rough boundary. These numerical results are very significant in understanding the impact of complex boundary conditions on hydrodynamic properties in nanofluidic theory and the design of nano-devices.

Methane has emerged as a key source of energy supplement and a vital greenhouse gas 1 . Understanding the hydrodynamic properties of methane is of significant interests from both theoretical research and nanofluidic device perspectives. Over the last several decades, numerical modeling of fluid flow inside micro and nanochannels with different boundary conditions plays a significant role in the optimal design of micro and nanofluidic devices, such as shale gas storage 2 , water purification 3 , drug delivery 4 , and nanomanufacturing 5 . The natures of hydrodynamic properties involved into these devices are predominated by the complex wall-fluid interaction force because of the microstructure of wall surfaces at the micro and nanoscales. Moreover, the continuum fluid theory is incompletely satisfied in nanoscale confined problems, where quantities such as velocity profiles do not remain classical Navier-Stokes (N-S) hydrodynamics [6][7][8] . Both numerical simulation and experimental studies indicated that the classical hydrodynamics theory was not valid when downsizing to the four times diameters of argon atom 7 or methane molecule [9][10][11] . Meanwhile, taking into consideration realistic interest of methane due to the fact that it is an important source of energy and a vital greenhouse gas 12 . Besides, to the best of our knowledge, investigation of the impact of nonslip and rough boundaries on the hydrodynamic properties of nanofluidic is rare. In this study, therefore, we propose the mixture boundary Poiseuille flow model to investigate the influence of boundary conditions on the hydrodynamic properties and explore the molecular mechanism of the confined methane nanofluidic.
Previous investigations of the hydrodynamic properties for nanofluidic were reported by Non-equilibrium molecular dynamic (NEMD) or equilibrium molecular dynamic (MD) simulations for various confined nanofluidics with different atom walls. The momentum transport characteristics of fluid molecules were reviewed by Cao et al. 6 for micro and nanofluidics at various nanochannel surfaces in micro-/nano-electro-mechanical systems. Some computational simulation techniques were discussed by Xie et al. 13 to study new thermal physical transport phenomena in the length range from nanoscale to micron. The origin of slip or nonslip boundary conditions was reported by MD simulations focusing on the electrostatic systems, pointing that the electrostatic interaction forces play a significant role for slip or nonslip boundary of the clay nanometer surfaces 14 . Meantime, in order to capture the microscopic details of the heat conduction accurately, Li et al. 15 developed a hybrid Monte Carlo computational framework to investigate large scale heat conduction mechanisms for ballistic-diffusive systems. Hydrodynamic and structural properties of the argon nanofluidic were systematically studied by Sofos et al. [16][17][18] using NEMD simulations. They showed that the transport behavior 16,17 and velocity profiles 18 presented asymmetric to the channel centerline within a grooved and a ribbed wall channel. Aminfar et al. 19,20 employed MD simulation to investigate the atomic behavior and nanoparticles aggregation of liquid-solid nanofluidic flows inside nanochannels, and pointed out that the microstructure and aggregation of nanoparticles were significantly influenced by the interaction force of particles at a specified body driving force 20 . Toghraie et al. 21 carried out MD simulation to investigate the agglutination characteristic of nanoparticles in a nanochannel. They showed that the agglutination time of Copper nanoparticles is faster than that of Platinum nanoparticles in nanochannel. Subsequently, the group of Toghraie [22][23][24][25] studied the effects of geometrical parameters, various nanochannel wall temperatures and the number of nanoparticles on the local physical property evolution and diffusion characteristic of the confined nanifluidic by using MD simulations. They found that the rough boundary conditions significantly influenced on the molecular mechanism and flow characteristic of Poiseuille flow in a nanochannels. He et al. 26 used the MD simulation to investigate the mass flux of methane nanofluidic in a rough nanoporous under the pressure gradient. Zhang [27][28][29] presented the flow factor approach model to investigate the hydrodynamics properties for nanometer fluid confined in different wall widths and various wall-fluid interactions, and showed that the interaction forces of wall-fluid influenced the hydrodynamics characteristics of confined nanofluidic significantly 28 . The dissipative particle dynamics (DPD) framework was applied to explore the impact of wall cavitations shape on flow feature of micron scale confined fluid by Kasiteropoulou et al. 30 . The study indicated that the density, temperature and pressure approached almost a constant in the center domain of nanochannel and their trait nearby the nanochannel surfaces relied on the cavitations size (or roughness). The influences of mass transport were studied by Yan et al. 31 using a Lattice Boltzmann method for volatile organic confined inside rough micropore, and indicated that the structure of micropore surfaces had significantly impact on the characteristics of mass transfer. Yu et al. 32 developed a multiscale Lattice Boltzmann method to explore the transport features for shale gas confined inside rough nanopores. However, the nature characteristics of wall materials were ignored in these researches. To resolve the complex wall-fluid interaction problem, NEMSMD framework was developed by Jiang et al. 33 to probe into the transport properties, the structural characteristics and flow behaviors for the methane fluid inside the various solid atom wall surfaces. They found that the nanochannel surfaces and interaction force between wall atoms and fluid molecules played a very important role on the nanofluidic flow 9,11,33 , and indicated that the parabolic characteristics between diffusion coefficient values and the inverse of nanochannel width were presented for methane fluid confined inside rough nanometer wall surfaces 9 .
Besides, many research works were implemented for investigating the slip length of methane (shale or natural gas) particularly confined in different nanopores, such as physical experiments 34 and numerical simulations 1,32,33,[35][36][37][38] . The slip velocities of methane flow confined in nanochannel were investigated by MD simulations, and showed that the influence of rough nanometer pore surfaces on the slip velocity was rather remarkable in an extremely small Kerogen nanopore 36,39 . The permeability and diffusivity of shale (methane and water) were explored by Chen et al. 40 using the lattice Boltzmann method for the reconstructed shale, and indicated that the Knudsen diffusion influenced significantly the transport mechanisms of shale gas through the porous geometrical configurations. The transport characteristics of methane inside nanofluidic were studied by Nan et al. 35 using the NEMD simulation, showed that the nanochannel size played a vital role in the slip length of confined nanofluidic when the pressure decreased to 10 MPa. Mirzaeifard et al 38 used multiscale modeling to study the influence of interface on the hydrodynamics characteristics for water and methane mixture systems, demonstrating that the interfacial tension decreased slightly with pressure drop or temperature increases when the fluid molecules approached the interface. Understanding the hydrodynamics characteristics has significant interests from both numerical simulation technique and fundamental perspectives for shale gas (Shale gas was a borne remarkable resource 1 ). With the develop of numerical simulation technique, methane adsorption behavior confined by organic nanochannel was investigated by Cao et al. 41 using Monte Carlo and MD technique. It has been presented that the surface adsorption plays a key contribution on the hydrodynamics features of nanofluidic flow inside nanopore 37 . Besides, the slip and nonslip boundary has been generally employed to explore flow characteristic of conventional methane nanometer confined fluid reservoirs 36 . The methane clathrate hydrate and adsorption in nanoporous materials are two independent methods for the methane storage technique in high density and pressure 42 . Majumder et al. 43 proposed that the slippage of nanopore wall surfaces had an outstanding influence on the flow behavior in carbon nanotubes. For nature methane, it would be of remarkable theoretical interest and realistic value to investigate hydrodynamic characteristics through the nature physical and asymmetric nanochannel. For example, its upper surface comprises nonslip boundary, and the lower surface adopts the nature physical boundary.
Investigating the effects of nonslip and rough boundary on the hydrodynamic behavior of confined methane nanofluidic, one of the new contributions of this study, which may be said to have not been reported in any of the previous work. This is one of the most important issues that the previous researchers have simply eliminated. Moreover, direct experiment investigation of the effect of the mixture boundary (the upper plate is nature physical

Modeling and simulation detail
Mixture boundary Poiseuille flow modeling. In this study, the mixture boundary Poiseuille flow modeling is conducted for dense confined methane nanofluidic using NEMSMD simulation (Fig. 1a). In order to investigate the impact of the nonslip boundary and rough boundary on hydrodynamic properties of confined methane nanofluidic in one simulation, the proposed model coupled the advantage of periodic Poiseulle flow model 44   www.nature.com/scientificreports/ molecule crosses the centerline (the details see Fig. 1a), the stress and velocity profiles for methane fluid with symmetrical boundary conditions are given by 44,45 In the above equations, τ zx , f ex , ρ , v X , η and D stand for the shear stress, body driving force, flow density, velocity in streaming, dynamic viscosity and the half of distance between the rough silicon atom plates, respectively. Meantime, the ABAB stacking is applied to construct the silicon atom for nanochannel plates due to the fact that it is an equilibrium configuration 46 . To quantify the roughness of natural physical wall surfaces, the rough nanochannel surfaces are constructed by using the periodic rectangular wave 33 where A is the amplitude, indicates the wavelength.
To explore the effect mechanism of boundary condition on the methane nanofluidic, the nonslip boundary is applied in the periodic Poiseuille flow, and the distance between two rough plates is 2D = 30.93 σ (σ = 4.01 Å) and the rough configuration ( = 4.31 σ , A=1.07 σ ) is adopted. In driving force direction (i.e. x-direction), an integer times of wavelength is employed to match the boundary condition of its periodicity. All interaction forces between methane molecules are conducted by our previous modification of optimized potential for liquid simulation (MOPLS) model in this study, because it could estimate accurately the mass transfer behaviors and structural characteristic for methane nanofluidic by the previous study of Jiang et al. 47 . It is defined as where ε and σ denote the energy well of methane molecule and length unit, respectively, q is point charge, ε 0 indicates the permittivity of vacuum, and r ij denotes the interaction distance between atoms (C or H) of methane, a and b are employed to distinguish the atom of C and H, α and β mark different methane molecules. In methane model: l CH = 1.087 Å and q C = −4q H = −0.572 e(e = 4.803 × 10 −10 esu ). The details of MOPLS model are listed below (see Table 1).
The NEMSMD framework is employed to manipulate the interaction force between wall atom and methane molecule (refer to Fig. 1), i.e. it is determined by coupling the CG potential of methane and silicon atom potential based on L-B mixing rule when the methane molecule arrives near the wall ( ≤ r cut ). The mainly advantage not only observe atomic information using the NEMSMD simulation, but also can ensure the interaction between wall atom and fluid molecule approach the realistic interaction of wall-fluid. The validity of NEMSMD framework can be referred to literature 10,33 . Furthermore, the CG potential model should be consistent with the potential function of silicon. In this paper, the LJ (12-6) potential is applied to calculate the interaction between two silicon atoms. Thus, the interaction potential function between methane molecule and wall atom is showed by where the parameter ε MS and σ MS are calculated by L-B mixing rule In this equation, the interaction factor χ = 1.5 is employed in all simulation cases. The interaction potential parameters between two silicon atoms are ε S = 1.6885 kJ/mol and σ S = 3.826 Å , the corresponding details can be seen in the literature 46 . The parameters of CG methane potential σ CG and ε CG are optimized by the relative entropy minimization framework 48 for all atom MOPLS in the bulk ensemble of methane under the selected state conditions. The details of optimized parameter for methane CG potential can be referred to our previous literature 49 . The parameter details are showed in Table 2. www.nature.com/scientificreports/ Simulation details. An equilibrium configuration ABAB stacking is employed for rough silicon atom walls, all atoms are conducted by the harmonic potential site r eq 50 , here r(t) is the position coordinate of silicon atom at time t . The second derivative value of silicon atom potential at r = r 0 = 2 1 / 6 σ S is used to determine the k w = 72ε S / 2 1 / 3 σ 2 S . In the direction of x coordinate axis, the body driving force ( f ex = 0.075 (ε/σ ) ) of methane molecules is employed to refrain the compressibility effects within the linear regime for confined methane nanofluidic 10,51 . The cutoff radius ( r cut = 2.5 σ ) is used to calculate the interaction force of different particles for all simulations. Moreover, the fourth Predictor Corrector method is employed to manipulate the motion of mass center for methane fluid molecule. The NVT ensemble is used to command the temperature of simulation system by Nosé-Hoover thermostat 45 . In order to control the simulations system's temperature and not effect on the flow dynamics, the thermostat is coupled to the thermostat in y and z directions, i.e. mv 2 y 2 = k B T 2 and mv 2 z 2 = k B T 2 , which can be referred to Bhadauria's research 52 (where m , v y , v z , k B and T are the methane molecule mass, velocity of y and z directions, Boltzmann constant and simulation temperature, respectively). Although it is not coupled to the velocity of x direction explicitly, this technique can ensure that the thermostat effect is obtained in each direction through intermolecular interactions. In this study, σ = 4.01 Å and ε k B = 142.87 K are the reduced unit, respectively. The total calculation steps are 8.0 × 10 5 . The time step is t = 0.001 ( 1.3 × 10 −15 s ). To collect the flow velocity, atom number density and temperature distributions of the confined methane nanofluidic accurately, the original calculation steps 3.0 × 10 5 are discarded. Furthermore, the motion information of methane fluid molecule is sampled by next 5.0 × 10 5 calculation steps. The velocity profiles and atom (C and H) number density distribution are collected via dividing the range in z-direction ( 2D + 2A ) by 900 bins. The velocity profiles are collected by 10,16 in order to calculate the velocity profiles of the nonslip boundary ( v N,x (z 1 ) ) and rough boundary ( v R,x (z 2 ) ) by using results of Eq. (8), the coordinate transformation methods on the local symmetry is given by The number density profiles of C and H atoms are determined by 33 and the temperature profiles across the fluid regions are calculated by 33 where v and v i indicate the velocity of macroscopic flow and velocity of molecule i, respectively. N is the total molecule number.
The shear viscosity is determined by the Poiseuille flow method for the methane nanofluidic, Eq. (2) can be reformulated as follows 10 www.nature.com/scientificreports/ in addition, the streaming velocity profile v x (z) is rewritten as the local viscosity is calculated by and the shear viscosities for the whole fluid domain of nonslip or rough boundary are determined by in the above Eqs. (13)(14)(15)(16)(17), η denotes the shear viscosity, v 0 and z 0 represent the maximum velocity and its position, and k is calculated via the fitted local velocity information derived from the proposed model. It is noted that the fitted local velocity technique will breakdown while the fluid molecules are confined inside the small width (less than 4 times of molecule diameters, the classical N-S equation fails for simple fluids 7,9 ) of nanochannels. In this study, we divided 26 layers in the z-direction (Fig. 1a) to calculate the shear viscosity and corresponding local values for methane fluid. Moreover, the collected local velocity samples are shown in Fig. 1b for calculating the local viscosity of the nonslip ( Num = 12 ) and rough ( Num = 14 ) boundary conditions in detail, i.e., the local viscosity of the conducted domain of the nonslip boundary is determined by the fitted local velocity (see from the top of Fig. 1b), and local shear viscosity of the conducted domain of the rough boundary is calculated via the fitted local velocity profiles (refer to the bottom of Fig. 1b). Finally, to study molecular mechanism of confined methane nanofluidic flowing the nonslip and natural physical (silicon atomic plates) boundary conditions, the whole nanochannel is divided into n x × n y × n z = (36 × 36 × 90) bins with each volume V bin = L x n x × L y n y × (L z + 4A) n z . The velocity field is calculated by 10 The three-dimensional distributions of translational and rotational kinetic energies are respectively calculated by the formula: and where m bin = mn bin , I x , v bin and w x, bin are mass, one component of the (diagonal) inertia tensor ( x indicate the number for x, y and z coordinate directions), average velocity, and average angular velocity with corresponding coordinate directions. And v i denotes the i-th molecule's velocity, n bin is the number of molecule in the bin. All the codes of program are developed by C++ based on Literature 45 , running on windows operating systems.

Results and discussion
To verify and examine the validity of mixture boundary Poiseuille flow model, the width of nanochannel is selected properly. Meanwhile, the atom number density, velocity and shear stress profiles are calculated. Moreover, the slip length is calculated near the rough boundary. To explore the impact mechanism of the asymmetric boundary condition (with nonslip and rough boundary) on methane nanofluidic, the local shear viscosity, the translational and rotational energies are discussed.
Density distribution, shear stress, velocity profile and slip length. In order to show the validity and availability of the proposed modeling in studying the hydrodynamic properties for confined methane nanofluidic, the width of nanochannel is selected properly 7,55 . Figure η lay = −ρf ex 2k lay ,   Figure 2d3 presents the velocity profiles agree with the analytical solution of the Poiseuille flow, which demonstrate that the distance of between two rough nanochannel walls is larger than 12 molecules diameter (about 4.968 nm ) if you want to calculate the viscosity by the fitted velocity profile method using the proposed modeling. Therefore, in the study, the distance 2D = 30.93 σ (about 12 nm ) between two rough nanochannel walls is selected to study the impact of boundary on the hydrodynamic properties of methane nanofluidic. The C and H atoms number density distribution profiles of methane are showed in Fig. 3 by NEMSMD simulations for different state points. It is obvious that the rough nanochannel surfaces and nonslip boundary have a significant influence on C and H atoms number density distribution when the fluid molecule moves to the boundary domain (nonslip and rough boundary). For different state points, the atoms number density distributions present strong oscillations near the rough nanochannel surfaces, and the methane fluid molecules remain extended period in some layers parallel to the nanochannel surfaces. These layers of the atom number density profiles are situated at the peaks visible because the methane fluid molecule suffered from the stronger repulsion force near the rough nanochannel surfaces ( ≤ r cut ). Moreover, the local atom (C or H) number density reduces with the temperature increase (or the density decrease) while they go into the fluid domain near the nonslip boundary. The values of atom number density approach a certain constant in the center of fluid domain between nonslip and rough boundary. Comparison of the H atom number density with that of C atom is showed. The positions of first peak value for the H and C atom number distribution indicate that the H atom can be firstly www.nature.com/scientificreports/ observed in the domain of near the nanochannel surfaces. The positions of the second valley and peak for the C atom and that of the third peak and valley for H atom present something similarly. Furthermore, the effect of the nonslip boundary on C atom number density is not more obvious than that of H atom, because the H atom is four times as many as C atom. This result also indicates that the spatial structure of methane is central symmetry spatial tetrahedral. Besides, the number density profile of atom (H or C) near the nonslip boundary is smaller than that of center fluid regions. The reason may be that the methane fluid molecule suffers from different shear forces (or shear-thinning). And the significant degree of this phenomenon depends on the state points, i.e., it increases with density decreasing (or temperature increasing). The shear stress profiles (reduced unit ε σ 3 ) are showed in Fig. 4 for selected state points based on the mixture boundary Poiseuille flow model using the NEMSMD simulations. Because of the powerful repulsion interaction force between the fluid molecule and silicon atom, the shear stress presents significant vibration for all cases in the fluid region near the nanochannel surfaces. Furthermore, with the temperature increasing, the shear stress slightly increases as the fluid molecules approach the domain of the nonslip boundary. The reason may be that the methane molecules are conducted by the same value of body driving force with opposite direction while it crosses the nonslip boundary. Besides, the thermal motion of methane molecules enhances with the system temperature increasing. However, we find that the shear stress is very close to the theoretical values inside the fluid domain between nonslip boundary and rough nanochannel wall. The similar result can be referred to the research of Backer et al. 44 for LJ fluid.
The temperature profiles of methane nanofluidic are calculated by the mixture boundary Poiseuille flow model using the NEMSMD simulations and plotted in Fig. 5 for different state points by Eq. (12). From the Fig. 5, the local temperature of the confined methane nanofluidic are not uniform everywhere, i.e., it shows that the fluid temperature near the nonslip boundary is higher than that of whole fluid regions. Moreover, the local temperature increases with the decreasing of the density (or increasing of the system temperature) when the fluid molecules arrive at the nonslip boundary. This reason may be that the fluid molecules are subject to the same body driving force with opposite direction while it crosses the nonslip boundary, which produces the latent heat of methane fluid molecule. Furthermore, the temperature profiles present much variation near the rough nanochannel surfaces, which are very close to the results of the DPD simulation by Kasiteropoulou et al. 30 for the Poiseuille flow confined inside grooved nanochannel. We ascribe this reason to the intricate interaction force (between wall atom and fluid molecule) which remarkably impacts the distributions of atom number density (Fig. 3). In addition, the values of temperature also illustrated a slight oscillation around the temperature of system in the www.nature.com/scientificreports/ center regions between nonslip boundary and rough nanochannel surfaces. The reason may be ascribed to the matter of fact that the temperature of the simulation system is regulated by calculating the kinetic energy in y and z coordinate directions. Besides, the statistical errors are not be ignored completely. These numerical results are generally similar to those of confined fluids by the NEMD simulations 53 . All numerical simulation results for temperature profiles of methane confined fluid also demonstrate that the mixture boundary Poiseuille flow model is reasonable and credible to study the confined methane nanofluidic using NEMSMD simulations. In order to illustrate the influence of the complex boundary on the hydrodynamic properties for the confined methane nanofluidic, the velocity profiles of the nonslip boundary are compared with that of rough boundary based on the local symmetrical and coordinate transformation methods. Figure 6a plots the velocity profiles of methane confined nanofluidic at the selected state points. As can be seen from the Fig. 6a, the irregular fluctuation of velocity values can be clearly observed near the rough nanochannel surfaces, this reason can be ascribed to the fact that the fluid molecules are subjected to the complex interaction force between wall atom and methane molecule inside the cavities of wall. Moreover, the velocity profiles of methane nanofluidic increase with the temperature increasing (or density decreasing), the velocity profiles demonstrate the central symmetry about the point of intersection for the z-position line and nonslip boundary. However, the velocity profiles of I(v x (z)) (III(v x (z)) ) and II(v x (z))(IV(v x (z)) ) are asymmetric. We ascribe it to the fluid molecule obtaining diverse   www.nature.com/scientificreports/ Shear viscosity. Table 3 shows the shear viscosity for the whole fluid conducted by nonslip and rough boundary conditions based on the proposed model with different state points (temperature and density). As it can be seen that the variation trend of shear viscosity is well displayed by the confined methane nanofluidic with the nonslip and rough boundary, the numerical results of the nonslip boundary are in better agreement with the experimental values than those of the fluid domain confined by rough boundary, which verifies the validity of the proposed model and the fitted local velocity technique in calculating shear viscosity for methane fluid. Their difference is attributed to different boundary conditions and the methane molecule model (to our knowledge, no molecule model can accurately predict all properties of methane at all state points 7,44,54 ). Moreover, the slip length of methane nanofluidic near the rough boundary may be dependent on the viscosity of fluid, i.e., the slip length decreases with the increasing of the viscosity (seen from Fig. 7 and Table 3). Last but not least, it is worthy of noting that the shear viscosity of the nonslip and rough boundary conditions are in agreement with the experiment values when the temperature is bigger than 170 K(or the density is less than 310.77 kg m 3 ). This fact is ascribed to the reason that the temperature (or the density) significantly influences the prediction results of shear viscosity for methane fluid. The numerical result also demonstrates the validity of fitted local velocity method in calculating the shear viscosity of confined nanofluidic. Furthermore, to better evaluate the influence of the boundary condition on the internal frictional force for confined methane nanofluidic, the local shear viscosity in different fluid layers is further depicted in Fig. 8. By comparing the local shear viscosity of the nonslip boundary with that of the rough boundary condition, it is observed that the present results regarding the local shear viscosity approach the experiment data for saturated liquid methane in the middle of fluid regions. The reason may be that the boundary conditions (nonslip boundary or rough boundary) obviously impact on the mobility of fluid molecule when the fluid molecule crosses the fluid region near the boundary, as are also shown by the atom number density distribution profiles in Fig. 3, stress profiles in Fig. 4 and velocity profiles in Fig. 6. In addition, the values of shear viscosity rapidly aggrandizes in the fluid nearby the rough nanochannel surfaces, which qualitatively agree with the simulation results for fluid confined in rough nanochannel walls 10,17 . And it differs significantly either in the nonslip boundary, as a matter of fact, the shear viscosity values slightly reduce in the fluid layer near the nonslip boundary, which is ascribed to the diminution of fluid layer molecule density (seen form Fig. 3

Molecular mechanism.
To obtain deeper understanding of the hydrodynamic properties dependence for the methane transport in nanochannel, we further investigate molecular mechanism of methane fluid flowing with nonslip and natural physical (silicon atomic plates) boundary conditions. However, there is no existing numerical and experimental data of the nanochannel confined flow for comparisons nonslip and rough boundary conditions in one simulation (the effect of the nonslip and rough boundary on the nanofluidic behaviors) so far. Figure 9 displays the streamline and dimensionless mainstream velocity (the unit of velocity is (ε/m) 1/2 ) contour of methane molecule using the proposed mixture boundary at the selected state points. It is observed that the streamlines are significantly distorted as the molecule is closed to the nonslip boundary and rough nanochannel surfaces. Actually, this flow behavior is that the methane molecule suffers from complex interaction force near the boundary regions. The fluid molecule is mainly subjected to wall atom interaction force in the fluid domain near the rough nanochannel surfaces. In the fluid domain of the nonslip boundary, however, the molecular movement of methane are significantly manipulated by the opposite driving force and intermolecular forces. As one can see from the bottom of Fig. 9 ( n z = 13, n z = 44 ), the disorder level of streamlines increases with the decreasing of temperature near the fluid domain for the rough and nonslip boundary. And the streamlines disorder level of near nonslip is more obvious than that of rough boundary. It indicates that the movement of fluid molecules near the boundary is not only manipulated by the velocity in x direction, but also the velocities in y and z direction should be considered. The impact of y and z directions' velocities on the movement fluid of molecules gradually increases with the decreasing of the distance between wall atom and fluid molecule. Besides, The impact of y and z directions' velocities on the movement fluid molecules nearby the rough boundary is larger than the nonslip boundary. These results demonstrate that the rough boundary has more significant influence on the molecular movement for confined nanofluidic than nonslip boundary. The above results of streamlines present that the velocity field suffer from the influence of boundary significantly, because of the complex interaction force near the boundary. Furthermore, as the temperature rises (or density decreases), the absolute value of mainstream velocity increases in the fluid region between nonslip boundary and rough nanochannel surfaces. It is attributed to two factors: the thermal motion of methane molecule increases when the temperature rises (or density decreases), and the viscosity of fluid decreases with the increasing of the temperature (decreasing of the density).
In addition, the above numerical results are mainly attributed to the rotational movement and translational mobility of methane molecules confined by the nonslip and rough boundary conditions. To verify and examine this point, the contour distributions of rotational and translational energies calculated via ensemble average are plotted in Figs. 10 and 11 based on the proposed model. As seen in Figs. 10 and 11, the smallest values of rotational and translational energies ( n z = 12 ) are observed near the rough nanochannel surfaces. The result is attributed to the reason that the fluid molecule obtained stronger constraining force from wall atom when it flow the domain near the nanochannel surfaces. And this constraining force slows down the translational motion and rotary movement of methane molecule. Actually, in the fluid region of near nanochannel surfaces ( ≤ r cut ), the interaction force of fluid molecules is conducted by MOPLS model ( U MOPLS r ij ), and that between wall atoms and fluid molecules is determined using the wall-fluid interaction potential ( U MS r ij ). However, Fig. 10 shows that the maximum values of rotational energy obviously observed near the nonslip boundary ( n z = 45 ) at state point ( T = 170 K and ρ = 310.77 kg m −3 ). This result is attributed to the moving in opposite directions of fluid molecules in the two sides of the nonslip boundary, leading to the thermal motion of molecules increasing. . We attribute to the rotational movement of fluid molecules stronger dependent on the state point, i.e., the higher the temperature is, the severer the rotational movement is. Furthermore, Fig. 11 presents that the translational kinetic energy shows the symmetry with nonslip boundary from the whole simulation domain, indicating the absolute value of velocity is almost symmetry with nonslip boundary. The contour of translational kinetic energy illustrates the asymmetric in the fluid region between nonslip boundary and rough boundary, which demonstrate that the different boundary conditions obviously influence on the translational kinetic energy of fluid molecule. Besides, the local rotational energy and translational kinetic energy www.nature.com/scientificreports/ obviously differ in the region near the nonslip and rough boundary conditions. The mainly reason ascribed to the fluid molecule suffered different interaction forces when it moves to the fluid regions near the nonslip and nature rough boundary. It is worthy of noting that the translational kinetic energy arrives at minimal value in fluid region near the nonslip boundary (the picture is also showed in Fig. 10). This result also indicates that the fluid molecules have stronger thermal motion near the nonslip boundary region than other fluid region. In nanoscale, the complex interaction force is the primary cause impacting on the hydrodynamic properties of confined fluid molecule. These numerical results demonstrate that the hydrodynamic characteristics of confined fluid strictly relay on the movement of fluid molecule, including the rotational motion and translational motion. In fact, the hydrodynamic characteristics of the confined fluid are strictly related to the state point of fluid and boundary condition. Therefore, the proposed model is credible to study the hydrodynamic properties of methane confined nanofluidic by the NEMSMD simulation in the extreme boundary condition (nonslip and rough boundary conditions).

Conclusions
The mixture boundary Poiseuille flow model is proposed to explore the influence of the rough and the nonslip boundary on hydrodynamic properties for methane confined nanofluidic in one NEMSMD simulation. All numerical results demonstrate that our proposed model is effective and credible scheme to study the extremely boundary confined fluid. The most major findings are obtained from this study also including: (1) The distance of between two rough nanochannel walls is larger than 12 molecules diameter (about 4.968 nm ) if you want to calculate the viscosity by the fitted velocity profile method using the proposed modeling.  www.nature.com/scientificreports/ (6) In the fluid domain near the boundary, the streamlines are distorted. Moreover, the rough boundary has more significant influence on the molecular movement for confined nanofluidic than nonslip boundary.
This study will provide a novelty understanding for the hydrodynamic behavior of methane nanofluidic with the nonslip and rough boundary conditions, and contribute to further theoretical research for methane nanofluidic evaluation and exploitation. The extension of this study for NEMSMD simulation according our previous works [9][10][11]22 affords researchers and engineers a good choice for simulations. Next work, the proposed model will be used to study the shale gas, water and methane mixture system, haze aerosol, the microelectronic design of nano-devices, and so on.

Data availability
All data generated or analyzed during this study are included in this published article.